home *** CD-ROM | disk | FTP | other *** search
SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) NNNNAAAAMMMMEEEE SSSSRRRROOOOTTTTMMMMGGGG, DDDDRRRROOOOTTTTMMMMGGGG - Constructs a modified Givens plane rotation SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS Single precision Fortran: CCCCAAAALLLLLLLL SSSSRRRROOOOTTTTMMMMGGGG ((((_d_1,,,, _d_2,,,, _b_1,,,, _b_2,,,, _r_p_a_r_a_m)))) C/C++: ####iiiinnnncccclllluuuuddddeeee <<<<ssssccccssssllll____bbbbllllaaaassss....hhhh>>>> vvvvooooiiiidddd ssssrrrroooottttmmmmgggg ((((ffffllllooooaaaatttt *_d_1,,,, ffffllllooooaaaatttt *_d_2,,,, ffffllllooooaaaatttt *_b_1,,,, ffffllllooooaaaatttt _b_2,,,, ffffllllooooaaaatttt *_r_p_a_r_a_m))));;;; Double precision Fortran: CCCCAAAALLLLLLLL DDDDRRRROOOOTTTTMMMMGGGG ((((_d_1,,,, _d_2,,,, _b_1,,,, _b_2,,,, _r_p_a_r_a_m)))) C/C++: ####iiiinnnncccclllluuuuddddeeee <<<<ssssccccssssllll____bbbbllllaaaassss....hhhh>>>> vvvvooooiiiidddd ddddrrrroooottttmmmmgggg ((((ddddoooouuuubbbblllleeee *_d_1,,,, ddddoooouuuubbbblllleeee *_d_2,,,, ddddoooouuuubbbblllleeee *_b_1,,,, ddddoooouuuubbbblllleeee _b_2,,,, ddddoooouuuubbbblllleeee *_r_p_a_r_a_m))));;;; IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN These routines are part of the SCSL Scientific Library and can be loaded using either the ----llllssssccccssss or the ----llllssssccccssss____mmmmpppp option. The ----llllssssccccssss____mmmmpppp option directs the linker to use the multi-processor version of the library. When linking to SCSL with ----llllssssccccssss or ----llllssssccccssss____mmmmpppp, the default integer size is 4 bytes (32 bits). Another version of SCSL is available in which integers are 8 bytes (64 bits). This version allows the user access to larger memory sizes and helps when porting legacy Cray codes. It can be loaded by using the ----llllssssccccssss____iiii8888 option or the ----llllssssccccssss____iiii8888____mmmmpppp option. A program may use only one of the two versions; 4-byte integer and 8-byte integer library calls cannot be mixed. The C and C++ prototypes shown above are appropriate for the 4-byte integer version of SCSL. When using the 8-byte integer version, the variables of type iiiinnnntttt become lllloooonnnngggg lllloooonnnngggg and the <<<<ssssccccssssllll____bbbbllllaaaassss____iiii8888....hhhh>>>> header file should be included. DDDDEEEESSSSCCCCRRRRIIIIPPPPTTTTIIIIOOOONNNN These routines compute the elements of a modified Givens plane rotation matrix. See the NOTES section of this man page for information about the interpretation of the data types described in the following arguments. PPPPaaaaggggeeee 1111 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) These routines have the following arguments: _d_1 First diagonal element. (input and output) SSSSRRRROOOOTTTTMMMMGGGG: Single precision. DDDDRRRROOOOTTTTMMMMGGGG: Double precision. On input, this value is the first diagonal element of the scaling matrix _D. On the first call to the routine, this value is typically 1.0. Subsequent calls typically use the value from the previous call. On output, this value is the first diagonal element of the updated scaling matrix _D''''. For C/C++, a pointer to this value is passed. _d_2 Second diagonal element. (input and output) SSSSRRRROOOOTTTTMMMMGGGG: Single precision. DDDDRRRROOOOTTTTMMMMGGGG: Double precision. On input, this is the second diagonal element of the scaling matrix _D. On the second call to the routine, this value is typically 1.0. Subsequent calls typically use the value from the previous call. On output, this value is the second diagonal element of the updated scaling matrix _D''''. For C/C++, a pointer to this value is passed. _b_1 x-coordinate. (input and output) SSSSRRRROOOOTTTTMMMMGGGG: Single precision. DDDDRRRROOOOTTTTMMMMGGGG: Double precision. On input, this value is the _x-coordinate of the vector used to define the angle of rotation, before scaling (multiplying by the matrix _D). On output, this value is the _x-coordinate of the rotated vector, before scaling (multiplying by the matrix _D''''). For C/C++, a pointer to this value is passed. _b_2 y-coordinate. (input) SSSSRRRROOOOTTTTMMMMGGGG: Single precision. DDDDRRRROOOOTTTTMMMMGGGG: Double precision. On input, this value is the _y-coordinate of the vector used to define the angle of rotation, before scaling (multiplying by the matrix _D). It is unchanged on output. _r_p_a_r_a_m Real array of dimension 5. (output) SSSSRRRROOOOTTTTMMMMGGGG: Single precision array. DDDDRRRROOOOTTTTMMMMGGGG: Double precision array. This array contains rotation matrix information. The routine sets up the computed elements in _r_p_a_r_a_m from inputs _d_1, _d_2, b1, and _b_2. PPPPaaaaggggeeee 2222 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSttttaaaannnnddddaaaarrrrdddd GGGGiiiivvvveeeennnnssss RRRRoooottttaaaattttiiiioooonnnn A standard Givens rotation (see SSSSRRRROOOOTTTTGGGG(3S)) is based on an orthogonal matrix _G that rotates points on a Cartesian _x_y-coordinate plane. To calculate the rotation matrix, you must provide the angle of rotation desired, or, equivalently, a vector (point) that lies along the angle of rotation. For a given planar point (_x_r, _y_r), _G is formed so that: _ _ _ _ _ _ _ _ | x' | | c s | | x | G | x | | | | | | r | | r | | 0 | = |-s c | | y | = | y | | | | | | r | | r | - - - - - - - - where _x' = sqrt (_x_r + _y_r) . With this rotation matrix _G, you can then convert any number of existing planar points to the new (rotated) _x_y-coordinate system. For _n points, the rotations would be as follows: _ _ _ _ _ _ | x | | c s | | x | | i | | | | i | | 0 |<- |-s c | | y | | i | | | | i | - - - - - - for _i = 1, 2, ..., _n MMMMooooddddiiiiffffiiiieeeedddd GGGGiiiivvvveeeennnnssss RRRRoooottttaaaattttiiiioooonnnn The algorithm for these routines is based on the following observation. The rotation matrix _G can be factored into a _s_c_a_l_i_n_g matrix (diagonal matrix) and _m_o_d_i_f_i_e_d rotation matrix _H, for which either the diagonal or the off-diagonal elements are units (that is, +1 or -1). Thus, to perform _m modified (scaled) rotations on _n planar points, requires only 2_n_m, rather than 4_n_m multiplications for the standard rotation. Because you may want to perform several successive rotations, this routine assumes that you have leftover scaling factors from your previous modified Givens rotation; that is, the routine requires you to input not only a planar rotation vector (_b_1, _b_2) but also the squares of the diagonal elements of the scaling matrix, _d_1 and _d_2. The actual rotation vector is specified as follows: _ _ _ _ _ _ 1/2_ _ | xr | | sqrt(d1) 0 | | b1 | D | b1 | | yr | = | 0 sqrt(d2) | | b2 | = | b2 | - - - - - - - - PPPPaaaaggggeeee 3333 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) where _d_1 and _d_2 are the input scaling factors. Given these inputs, the routine generates a new modified rotation matrix _H with units for either the diagonal or off-diagonal elements, and new elements _d_1' and _d_2' for the new scaling factors, and rotated but unscaled vector (_b_1', 0), with the following results: _ _ 1/2_ _ 1/2 _ _ 1/2_ _ _ _ G | xr | G D | b1 | D' H | b1 | = D' | h11 h12 | | b1 | | yr | = | b2 | = | b2 | | h21 h22 | | b2 | - - - - - - - - - - 1/2_ _ _ _ D' | b1' | | x' | = | 0 | = | 0 | - - - - where: 1/2 D' ( = diag {sqrt(d1'), sqrt(d2')} ) uses the updated scaling factors: * _d_1'''' * _d_2'''', which are _d_1 and _d_2 on output. _H is stored in the output array argument _r_p_a_r_a_m; _b_1' is stored as _b_1 on output. _D''''1111////2222 HHHH equals _G _D _1/_2, nnnnooootttt _G, as implied earlier. You must account for the old scaling factors when calculating the new scaling factors. After calculating the matrix _H by using SSSSRRRROOOOTTTTMMMMGGGG/DDDDRRRROOOOTTTTMMMMGGGG, you can then use it in SSSSRRRROOOOTTTTMMMM/DDDDRRRROOOOTTTTMMMM to convert points to the new coordinate system. MMMMeeeeaaaannnniiiinnnngggg ooooffff tttthhhheeee OOOOuuuuttttppppuuuutttt VVVVaaaalllluuuueeeessss The output values are returned through arguments _d_1, _d_2, _b_1, and _r_p_a_r_a_m. The scaling factors _d_1 and _d_2 are updated with each call to the routine. Although SSSSRRRROOOOTTTTMMMM/DDDDRRRROOOOTTTTMMMM does not need the updated factors, they are needed in two other important contexts: * As input for subsequent calls to this routine. * As scaling factors for rotated but unscaled points (_x_i, _y_i), which are output from SSSSRRRROOOOTTTTMMMM./DDDDRRRROOOOTTTTMMMM. PPPPaaaaggggeeee 4444 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) In this second usage, the actual (scaled) points would be given by (sqrt(d1)*x(i), sqrt(d2)*y(i)) Doing this operation frequently on all your points is counterproductive. The main advantage of the modified rotation algorithm is to reduce the number of operations. If you fold in the scaling factors after each rotation, you are performing the same number of operations as in the standard Givens rotation. These two uses for the scaling factors are mutually exclusive; that is, if you fold the scaling factors back into all your points, you no longer need those factors for these routines. After folding in the scaling factors, and before the next call to SSSSRRRROOOOTTTTMMMMGGGG/DDDDRRRROOOOTTTTMMMMGGGG, reset _d_1 and _d_2 to 1.0. On output, _b_1 represents the new _x-coordinate after rotating (but before scaling) the rotation vector. Although the _y-coordinate of this vector is 0.0 (see the the previous discussion), the corresponding value _b_2 is unchanged on output. The output array argument _r_p_a_r_a_m specifies the format of matrix _H, and it holds the nonunit values of _H. This is the only output of these routines that SSSSRRRROOOOTTTTMMMM/DDDDRRRROOOOTTTTMMMM requires. Each element of _r_p_a_r_a_m has a specific meaning, as follows: _r_p_a_r_a_m(1) a flag parameter that specifies how the matrix is stored = 0.0 Off-diagonal elements of _H are units. = 1.0 Diagonal elements of _H are units. = -1.0 Rescaling case (see the following subsection). = -2.0 _H is the identity matrix; no rotation needed. _r_p_a_r_a_m(2) = _h_1,_1 if needed _r_p_a_r_a_m(3) = _h_2,_1 if needed _r_p_a_r_a_m(4) = _h_1,_2 if needed _r_p_a_r_a_m(5) = _h_2,_2 if needed CCCCaaaallllccccuuuullllaaaattttiiiinnnngggg tttthhhheeee OOOOuuuuttttppppuuuutttt VVVVaaaalllluuuueeeessss The following presents the algorithm for calculating the output values _d_1, _d_2, _b_1, and _r_p_a_r_a_m, based on the input values _d_1, _d_2, _b_1, and _b_2. This algorithm is presented without explanation or proof. For a more complete discussion of the modified Givens algorithm, see the papers listed in the SEE ALSO section. PPPPaaaaggggeeee 5555 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) CCCCaaaasssseeee 1111:::: bbbb2222 ==== 0000 ((((ttttrrrriiiivvvviiiiaaaallll ccccaaaasssseeee)))) In this case, no rotation is needed. The flag value, _r_p_a_r_a_m(1), is set to -2.0. When passed to routine SSSSRRRROOOOTTTTMMMM/DDDDRRRROOOOTTTTMMMM, this flag indicates that it should not do any rotations. On output from SSSSRRRROOOOTTTTMMMMGGGG/DDDDRRRROOOOTTTTMMMMGGGG, _d_1, _d_2, and _b_1 are unchanged. CCCCaaaasssseeee 2222:::: |||| ssssqqqqrrrrtttt((((dddd1111)))) **** bbbb1111 |||| >>>> |||| ssssqqqqrrrrtttt((((dddd2222))))****bbbb2222 |||| ((((||||xxxxrrrr|||| >>>> ||||yyyyrrrr||||)))) In this case, the diagonal elements of _H (h1,1 and h2,2) are set to 1.0. Thus, the _r_p_a_r_a_m values set on output are as follows: rparam(1) <-- 0.0 rparam(3) <-- h2,1 = -b2/b1 rparam(4) <-- h1,2 = (d2*b2)/(d1*b1) The output values of the scaling factors are as follows: d1 <-- d1' = d1/u d2 <-- d2' = d2/u where u = det(H) = 1 + (d2*b2*b2)/(d1*b1*b1) The output value of _b_1 is as follows: b1 <-- b1' = b1*u If rescaling is needed, SSSSRRRROOOOTTTTMMMMGGGG/DDDDRRRROOOOTTTTMMMMGGGG will further modify these output values before the end of the routine. See case 4 later in this subsection. CCCCaaaasssseeee 3333:::: |||| ssssqqqqrrrrtttt((((dddd1111))))****bbbb1111 |||| <<<<==== |||| ssssqqqqrrrrtttt((((dddd2222)))) **** bbbb2222 |||| ((((||||xxxxrrrr|||| <<<<==== ||||yyyyrrrr||||)))) In this case, the off-diagonal elements of _H are units (to be specific, _h_2,_1 = -1 and _h_1,_2 = 1). Thus, the _r_p_a_r_a_m values set on output are as follows: rparam(1) <-- 1.0 rparam(2) <-- h1,1 = (d1*b1)/(d2*b2) rparam(5) <-- h2,2 = b1/b2 The output values of the scaling factors are as follows: PPPPaaaaggggeeee 6666 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) d1 <-- d1' = d2/u d2 <-- d2' = d1/u where u = det(H) = 1 + (d1 * b1*b1)/(d2 * b2* b2) The output value of _b_1 is as follows: b1 <-- b1' = b2*u If rescaling is needed, SSSSRRRROOOOTTTTMMMMGGGG/DDDDRRRROOOOTTTTMMMMGGGG will further modify these output values before the end of the routine. See case 4. CCCCaaaasssseeee 4444:::: Rescaling If the scaling factors become either very large or very small, the scaling and rotation operations may lose a lot of accuracy; therefore, each scaling factor from case 2 or case 3 is kept within the range: gamma**(-2) <= | d(i)' | <= gamma**2 , for i = 1, 2 where d(1)' = d1', d(2)' = d2', and gamma = 2.0**12 = 4096.0. At the end of case 2 or case 3 assignments, if either of the scaling factors falls outside this range, this routine must rescale that factor (and the corresponding elements of _H) to bring its size back within the specified range (48 - log2(_g_a_m_m_a) = _4_8 - _1_2 = _3_6 _b_i_t_s (~ _1_0 _d_e_c_i_m_a_l _d_i_g_i_t_s)) Rescaling is performed as follows: If either _d_1 or _d_2 is 0, no rescaling is done; otherwise, let q(i) = int(log_base_gamma(sqrt(|d(i)'|))) = int(log2(|d(i)'|)/24) , for i = 1 , 2. Then the following is true: q(i) < 0 , if | d(i)'| < gamma**2 q(i) = 0 , if gamma**(-2) <= |d(i)'| <= gamma**2 q(i) > 0 , if | d(i)'| > gamma**2 PPPPaaaaggggeeee 7777 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) Furthermore, |q(i)| represents the number of times d(i)' must be multiplied (or divided) by _g_a_m_m_a**2 to return it to the proper range of values. In this case, the _r_p_a_r_a_m values set on output are as follows: rparam(1) <-- -1.0 rparam(2) <-- h1,1' = h1,1*gamma**q(1) rparam(3) <-- h2,1' = h2,1*gamma**q(2) rparam(4) <-- h1,2' = h1,2*gamma**q(1) rparam(5) <-- h2,2' = h2,2*gamma**q(2) The output values of the scaling factors are as follows: d1 <-- d1'' = d1'*gamma**(-2*q(1)) d2 <-- d2'' = d2'*gamma**(-2*q(2)) The output value of _b_1 is as follows: b1 <-- b1'' = b1'*gamma**q(1) NNNNOOOOTTTTEEEESSSS The following data types are described in this documentation: TTTTeeeerrrrmmmm UUUUsssseeeedddd DDDDaaaattttaaaa ttttyyyyppppeeee Fortran: Array dimensioned _n xxxx((((nnnn)))) Integer IIIINNNNTTTTEEEEGGGGEEEERRRR (IIIINNNNTTTTEEEEGGGGEEEERRRR****8888 for ----llllssssccccssss____iiii8888[[[[____mmmmpppp]]]]) Single precision RRRREEEEAAAALLLL Double precision DDDDOOOOUUUUBBBBLLLLEEEE PPPPRRRREEEECCCCIIIISSSSIIIIOOOONNNN C/C++: Array dimensioned _n xxxx[[[[_n]]]] Integer iiiinnnntttt (lllloooonnnngggg lllloooonnnngggg for ----llllssssccccssss____iiii8888[[[[____mmmmpppp]]]]) Single precision ffffllllooooaaaatttt Double precision ddddoooouuuubbbblllleeee PPPPaaaaggggeeee 8888 SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SSSSRRRROOOOTTTTMMMMGGGG((((3333SSSS)))) SEE ALSO IIIINNNNTTTTRRRROOOO____SSSSCCCCSSSSLLLL(3S), IIIINNNNTTTTRRRROOOO____BBBBLLLLAAAASSSS1111(3S) IIIINNNNTTTTRRRROOOO____CCCCBBBBLLLLAAAASSSS(3S) for information about using the C interface to Fortran 77 Basic Linear Algebra Subprograms (legacy BLAS) set forth by the Basic Linear Algebra Subprograms Technical Forum. SSSSRRRROOOOTTTT(3S), SSSSRRRROOOOTTTTGGGG(3S), SSSSRRRROOOOTTTTMMMM(3S) Gentleman, W. M., "Least Squares Computations by Givens Transformations Without Square Roots," _J_o_u_r_n_a_l _o_f _t_h_e _I_n_s_t_i_t_u_t_e _f_o_r _M_a_t_h_e_m_a_t_i_c_a_l _A_p_p_l_i_c_a_t_i_o_n_s 12 (1973), pp. 329 - 336. Lawson, C., Hanson, R., Kincaid, D., and Krogh, F., "Basic Linear Algebra Subprograms for Fortran Usage," _A_C_M _T_r_a_n_s_a_c_t_i_o_n_s _o_n _M_a_t_h_e_m_a_t_i_c_a_l _S_o_f_t_w_a_r_e, 5 (1979), pp. 308 - 325. PPPPaaaaggggeeee 9999 
